{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "54bd32c5",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-12-09T06:52:38.038063Z",
     "start_time": "2023-12-09T06:52:37.104554Z"
    },
    "execution": {
     "iopub.execute_input": "2024-02-19T12:37:25.770630Z",
     "iopub.status.busy": "2024-02-19T12:37:25.770630Z",
     "iopub.status.idle": "2024-02-19T12:37:26.547094Z",
     "shell.execute_reply": "2024-02-19T12:37:26.546589Z",
     "shell.execute_reply.started": "2024-02-19T12:37:25.770630Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.fft as fft\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from tqdm import tqdm\n",
    "import pandas as pd\n",
    "from cv2 import dct\n",
    "import imageio\n",
    "import quaternion as qt\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "57572914-634d-4f37-a95f-efb42cd8fff5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-19T13:18:22.035758Z",
     "iopub.status.busy": "2024-02-19T13:18:22.035758Z",
     "iopub.status.idle": "2024-02-19T13:18:22.078644Z",
     "shell.execute_reply": "2024-02-19T13:18:22.076746Z",
     "shell.execute_reply.started": "2024-02-19T13:18:22.035758Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def image_hash(img_path,k):\n",
    "    img = processing(img_path)\n",
    "    feature = image_feature(img)\n",
    "    h_i = gen_hashing(feature,k)\n",
    "    return h_i\n",
    "\n",
    "def processing(img_path):\n",
    "    \"\"\"\n",
    "    input：图片的路径\n",
    "    output：处理后的RGB图片\n",
    "    \"\"\"\n",
    "    try:\n",
    "        img = cv2.imread(img_path)\n",
    "        img = cv2.resize((img), (512,512), interpolation=cv2.INTER_LINEAR)\n",
    "    except:\n",
    "        img = imageio.mimread(img_path)\n",
    "        img = np.array(img)\n",
    "        img = img[0]\n",
    "        img = img[:, :, 0:3]\n",
    "        img = cv2.resize((img), (512,512), interpolation=cv2.INTER_LINEAR)\n",
    "#     out = cv2.GaussianBlur(img, (3, 3),1.3) # 使用python自带的高斯低通滤波\n",
    "    kernel = np.ones((7,7), np.float32) / 49 # 平均滤波器\n",
    "    out = cv2.filter2D(img, -1 , kernel=kernel)  # 二维滤波器\n",
    "    out = cv2.cvtColor(out, cv2.COLOR_BGR2RGB)\n",
    "    ring_img = ring_partition(out,256,256)\n",
    "    return ring_img.transpose(2,1,0)\n",
    "\n",
    "def ring_partition(img,n,image_block_size):\n",
    "    \"\"\"\n",
    "    使用环形分区形成二次图像\n",
    "    img:原始图片\n",
    "    n:环形分区\n",
    "    return:二次图像\n",
    "    \"\"\"\n",
    "    radius = [0]\n",
    "    u_a = image_block_size\n",
    "    radius.append(np.sqrt(u_a / np.pi))\n",
    "    R_result = np.zeros((n,u_a))\n",
    "    G_result = np.zeros((n,u_a))\n",
    "    B_result = np.zeros((n,u_a))\n",
    "    for i in range(1,n):\n",
    "        radius.append(np.sqrt((u_a + np.pi * radius[i] ** 2) / np.pi))\n",
    "    # radius = np.ceil(radius)\n",
    "    for i in range(len(radius)-1):\n",
    "        inner_radius = radius[i]\n",
    "        outer_radius = radius[i+1]\n",
    "        ring_result = ring_value(img,inner_radius,outer_radius,u_a)\n",
    "        R_result[i] = ring_result[0]\n",
    "        G_result[i] = ring_result[1]\n",
    "        B_result[i] = ring_result[2]\n",
    "    return np.array([R_result,G_result,B_result])\n",
    "\n",
    "def ring_value(img,inner_radius,outer_radius,u_a):\n",
    "    \"\"\"\n",
    "    获取环形分区中不为0的值，并且以升序的形式排序\n",
    "    return:(3,448)\n",
    "    \"\"\"\n",
    "    center = (img.shape[0],img.shape[0])\n",
    "    x,y = np.ogrid[0:img.shape[0],0:img.shape[0]]\n",
    "    distance = np.sqrt((x-img.shape[0]/2)**2 + (y-img.shape[0]/2)**2)\n",
    "    mask = (distance < outer_radius) & (distance >=inner_radius)\n",
    "    R = img[mask,0]\n",
    "    G = img[mask,1]\n",
    "    B = img[mask,2]\n",
    "    result = np.zeros((3,u_a))\n",
    "    for index,i in enumerate((R,G,B)):\n",
    "        xp = np.linspace(1, len(i), len(i))  \n",
    "        # 定义要插值的x值\n",
    "        x = np.linspace(1, len(i), u_a)\n",
    "        result[index] = np.sort(np.interp(x, xp, np.sort(i)))\n",
    "    return result\n",
    "\n",
    "def image_feature(img):\n",
    "    \"\"\"\n",
    "    iamge:(512,512,3)\n",
    "    return: array(512,512)\n",
    "    \"\"\"\n",
    "    rows,cols,_ = img.shape\n",
    "    quaternion_array = np.array([[np.quaternion(0, img[i,j,0],img[i,j,1],img[i,j,2]) for j in range(cols)] for i in range(rows)])\n",
    "    return np.fft.fftshift(fft2D(quaternion_array))\n",
    "    \n",
    "def exp_quaternion(u):\n",
    "    \"\"\"\n",
    "    四元数的指数运算\n",
    "    \"\"\"\n",
    "    # 手动计算四元数的模，并用它来规范化四元数\n",
    "    norm_u = np.abs(u)\n",
    "    v = u / (norm_u+(1e-80))\n",
    "    # 计算指数形式的指数运算\n",
    "    exponential_form = np.cos(norm_u) + v * np.sin(norm_u)\n",
    "    return exponential_form\n",
    "\n",
    "def fft1D(signal):\n",
    "    \"\"\"\n",
    "    signal:为一维的四元数序列\n",
    "    return:*\n",
    "    \"\"\"\n",
    "    # 定义单位纯四元数 u = ai + bj + ck，其中 |u| = 1\n",
    "    # 构建单位纯四元数\n",
    "    # 第一种四元数论文常用单位四元数\n",
    "    d1=3**(1/2)\n",
    "    u_lum = qt.quaternion(0, 1/d1, 1/d1, 1/d1)\n",
    "    # 第二种单位四元数\n",
    "    d2=68**(1/2)\n",
    "    u_perc=qt.quaternion(0, 0, -2/d2, 8/d2)\n",
    "    N = len(signal)\n",
    "    if N <= 1:\n",
    "        return signal\n",
    "    even = fft1D(signal[0::2])\n",
    "    odd = fft1D(signal[1::2])\n",
    "    #一维傅里叶变换公式\n",
    "    T = [exp_quaternion(-2*u_lum * np.pi * k / N) * odd[k] for k in range(N // 2)]\n",
    "    # T = [exp_quaternion(-2 * u_perc * np.pi * k / N) * odd[k] for k in range(N // 2)]\n",
    "    return np.concatenate([even + T, even - T])\n",
    "    \n",
    "def fft2D(image):\n",
    "    \"\"\"\n",
    "    image:四元数矩阵(256*256)\n",
    "    return:返回经过傅里叶变换后的四元数矩阵\n",
    "    \"\"\"\n",
    "    M, N = image.shape\n",
    "    if M <= 1 and N <= 1:\n",
    "        return image\n",
    "    # FFT along rows\n",
    "    rows = np.array([fft1D(row) for row in image])\n",
    "    # FFT along columns\n",
    "    cols = np.array([fft1D(col) for col in rows.T]).T\n",
    "    return cols\n",
    "\n",
    "def gen_hashing(feature_matrix,k):\n",
    "    \"\"\"\n",
    "    选取振幅系数矩阵中间偏左上的矩阵作为特征。\n",
    "    input:array (256,256)\n",
    "    output:list (x)\n",
    "    \"\"\"\n",
    "    coeff = []\n",
    "    start = int(feature_matrix.shape[0]/2-k//2+1)\n",
    "    for i in range(start,start+k):\n",
    "        for j in range(start,start+k):\n",
    "            coeff.append(feature_matrix[i][j].abs())\n",
    "    return np.array([1 if coeff[i] >= coeff[i+1] else 0 for i in range(len(coeff)-1)])\n",
    "    \n",
    "def dist_img(h1,h2):\n",
    "    return sum(np.abs(h1-h2))/len(h1)\n",
    "\n",
    "def dis_different_dir(path,k,des_path):\n",
    "    # 目录下图片之间的距离\n",
    "    if not os.path.exists(des_path):\n",
    "        os.makedirs(des_path)\n",
    "    dirs = os.listdir(path)\n",
    "    image_set_list = []\n",
    "    image_hashing_value_set = []\n",
    "    for i in tqdm(dirs,ncols = 50):\n",
    "    # for i in dirs:\n",
    "        image_hashing_value_set.append(image_hash(os.path.join(path, i),k))\n",
    "    for i in range(len(image_hashing_value_set)):\n",
    "        for j in range(i+1,len(image_hashing_value_set)):\n",
    "            image_set_list.append(dist_img(image_hashing_value_set[i],image_hashing_value_set[j]))\n",
    "    for i in range(int(np.ceil(len(image_set_list)/1000000))-1):\n",
    "        start = i*1000000\n",
    "        end = (i+1)*1000000\n",
    "        pd.DataFrame(image_set_list[start:end]).to_excel(os.path.join(des_path,f\"different_image_result_{i+1}.xlsx\"),index = False)\n",
    "    pd.DataFrame(image_set_list[end+1:]).to_excel(os.path.join(des_path,f\"different_image_result_{i+2}.xlsx\"),index = False)\n",
    "    return image_set_list,image_hashing_value_set\n",
    "\n",
    "def dis_similar_dir(path,k):\n",
    "    \"\"\"\n",
    "    path:相同图片每种操作的目录\n",
    "    k:选取中间多少元素作为哈希值\n",
    "    return:所有相同图片哈希距离的列表\n",
    "    \"\"\"\n",
    "    # 计算相同类型图片的距离\n",
    "    dir_rotation = os.listdir(path)\n",
    "    dis_similar = []\n",
    "    for i in dir_rotation:\n",
    "        if (os.path.splitext(i)[1] == \".bmp\") | (os.path.splitext(i)[1] == \".jpg\"):\n",
    "            h1 = image_hash(os.path.join(path, i),k)\n",
    "            dir_temp = os.path.join(path,os.path.splitext(i)[0])\n",
    "            for j in os.listdir(dir_temp):\n",
    "                h2 = image_hash(os.path.join(dir_temp,j),k)\n",
    "                dis_similar.append(dist_img(h1,h2))\n",
    "    return dis_similar\n",
    "\n",
    "def dir_dis_similar_dir(total_path,k,des_path):\n",
    "    \"\"\"\n",
    "    total_path:相同图片各种操作的根目录\n",
    "    des_path:存放结果的目录\n",
    "    return：存放每一种操作相同图片哈希距离的字典\n",
    "    \"\"\"\n",
    "    if not os.path.exists(des_path):\n",
    "        os.makedirs(des_path)\n",
    "    fold_each_attack = os.listdir(total_path)\n",
    "    dir_each_attack= {}\n",
    "    if not os.path.exists(os.path.join(des_path,\"similar_image_result\")):\n",
    "        os.makedirs(os.path.join(des_path,\"similar_image_result\"))\n",
    "    # total_path = os.path.join(total_path,\"similar_image_result\")\n",
    "    for i in tqdm(fold_each_attack,ncols = 50):\n",
    "        attack_path = os.path.join(total_path,i)\n",
    "        dir_each_attack[i] = dis_similar_dir(attack_path,k)\n",
    "        file_name = i + \".xlsx\"\n",
    "        pd.DataFrame(dir_each_attack[i]).to_excel(os.path.join(des_path,\"similar_image_result\",file_name),index = False)\n",
    "    return dir_each_attack\n",
    "\n",
    "def my_auc(image_similar, image_different):\n",
    "    from scipy import integrate\n",
    "    sum1 = len(image_similar)\n",
    "    sum2 = len(image_different)\n",
    "    tpr = []\n",
    "    fpr = []\n",
    "    threshold_max = (max(max(image_similar), max(image_different)))\n",
    "    for threshold in np.linspace(0,threshold_max,50):\n",
    "        tpr_number = len(image_similar[image_similar < threshold])\n",
    "        fpr_number = len(image_different[image_different < threshold])\n",
    "        tpr.append(tpr_number/sum1)\n",
    "        fpr.append(fpr_number/sum2)\n",
    "    return integrate.trapezoid(tpr,fpr)\n",
    "\n",
    "def stat_analysis(result_path):\n",
    "    \"\"\"\n",
    "    result_path:保存结果的路径\n",
    "    \"\"\"\n",
    "    dirs = os.listdir(result_path)\n",
    "    for i in dirs:\n",
    "        file_path = os.path.join(result_path, i)\n",
    "        if os.path.isdir (file_path):\n",
    "            for j in os.listdir(file_path):\n",
    "                xlsx_path = os.path.join(file_path,j)\n",
    "                data = pd.read_excel(xlsx_path,header = 0)[0]\n",
    "                print(f\"{j}的最小值：{min(data)}最大值：{max(data)}平均值：{np.mean(data)}\")\n",
    "        else:\n",
    "            xlsx_path = file_path\n",
    "            data = pd.read_excel(xlsx_path,header = 0)[0]\n",
    "            print(f\"{i}的最小值：{min(data)}最大值：{max(data)}平均值：{np.mean(data)}\")\n",
    "\n",
    "def auc_analysis(dir_path):\n",
    "    \"\"\"\n",
    "    dir_path:保存结果的路径\n",
    "    \"\"\"\n",
    "    similar_image_result_path = os.path.join(dir_path,\"similar_image_result\")\n",
    "    data_different =pd.Series()\n",
    "    for i in os.listdir(dir_path):\n",
    "        if not os.path.isdir(os.path.join(dir_path,i)):\n",
    "            print(os.path.join(dir_path,i))\n",
    "            temp = pd.read_excel(os.path.join(dir_path,i),header = 0)[0]\n",
    "            data_different = pd.concat((data_different,temp),axis=0)\n",
    "    similar_image_result_path_total = pd.Series()\n",
    "    for i in os.listdir(similar_image_result_path):\n",
    "        data_similar = pd.read_excel(os.path.join(similar_image_result_path,i),header = 0)[0]\n",
    "        similar_image_result_path_total= pd.concat((similar_image_result_path_total,data_similar),axis=0)\n",
    "        print(f\"{i}和不同图片的auc值：{my_auc(data_similar,data_different)}\")\n",
    "    print(f\"所有相同图片和不同图片 的auc值：{my_auc(similar_image_result_path_total,data_different)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "2db639b2-ef7f-4e81-83f0-32c3698e0e99",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-19T13:18:22.871332Z",
     "iopub.status.busy": "2024-02-19T13:18:22.870335Z",
     "iopub.status.idle": "2024-02-20T16:00:22.065914Z",
     "shell.execute_reply": "2024-02-20T16:00:22.063412Z",
     "shell.execute_reply.started": "2024-02-19T13:18:22.870335Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 9/9 [22:02:55<00:00, 8819.46s/it]\n",
      "100%|███████| 3307/3307 [4:32:19<00:00,  4.94s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "different_image_result_1.xlsx的最小值：0.09821428571428571最大值：0.6919642857142857平均值：0.370682205357143\n",
      "different_image_result_2.xlsx的最小值：0.08035714285714286最大值：0.6607142857142857平均值：0.36851799107142874\n",
      "different_image_result_3.xlsx的最小值：0.1116071428571429最大值：0.6696428571428571平均值：0.3711386964285713\n",
      "different_image_result_4.xlsx的最小值：0.1205357142857143最大值：0.6696428571428571平均值：0.36642249999999993\n",
      "different_image_result_5.xlsx的最小值：0.1026785714285714最大值：0.65625平均值：0.3667240089285715\n",
      "different_image_result_6.xlsx的最小值：0.09821428571428571最大值：0.6919642857142857平均值：0.3702814202566999\n",
      "brightness.xlsx的最小值：0.0最大值：0.1919642857142857平均值：0.055877976190476186\n",
      "contrast.xlsx的最小值：0.0最大值：0.1205357142857143平均值：0.033984375\n",
      "gamma correction.xlsx的最小值：0.0最大值：0.3169642857142857平均值：0.09026744378306878\n",
      "gaussian.xlsx的最小值：0.05357142857142857最大值：0.2946428571428572平均值：0.1806175595238095\n",
      "jpeg.xlsx的最小值：0.02232142857142857最大值：0.3258928571428572平均值：0.1557477678571429\n",
      "rotation.xlsx的最小值：0.0最大值：0.1383928571428572平均值：0.049641927083333336\n",
      "salt_pepper.xlsx的最小值：0.03125最大值：0.2946428571428572平均值：0.13908617424242425\n",
      "speckle.xlsx的最小值：0.05357142857142857最大值：0.2946428571428572平均值：0.1806175595238095\n",
      "watermark.xlsx的最小值：0.0最大值：0.01339285714285714平均值：0.0008184523809523809\n",
      "./新数据集结果/纯四元数版本_环形分区结果/different_image_result_1.xlsx\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\10047\\AppData\\Local\\Temp\\ipykernel_16736\\3643798627.py:245: FutureWarning: The behavior of array concatenation with empty entries is deprecated. In a future version, this will no longer exclude empty items when determining the result dtype. To retain the old behavior, exclude the empty entries before the concat operation.\n",
      "  data_different = pd.concat((data_different,temp),axis=0)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "./新数据集结果/纯四元数版本_环形分区结果/different_image_result_2.xlsx\n",
      "./新数据集结果/纯四元数版本_环形分区结果/different_image_result_3.xlsx\n",
      "./新数据集结果/纯四元数版本_环形分区结果/different_image_result_4.xlsx\n",
      "./新数据集结果/纯四元数版本_环形分区结果/different_image_result_5.xlsx\n",
      "./新数据集结果/纯四元数版本_环形分区结果/different_image_result_6.xlsx\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\10047\\AppData\\Local\\Temp\\ipykernel_16736\\3643798627.py:249: FutureWarning: The behavior of array concatenation with empty entries is deprecated. In a future version, this will no longer exclude empty items when determining the result dtype. To retain the old behavior, exclude the empty entries before the concat operation.\n",
      "  similar_image_result_path_total= pd.concat((similar_image_result_path_total,data_similar),axis=0)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "brightness.xlsx和不同图片的auc值：0.9999975505978569\n",
      "contrast.xlsx和不同图片的auc值：0.9999996070742788\n",
      "gamma correction.xlsx和不同图片的auc值：0.998894684363306\n",
      "gaussian.xlsx和不同图片的auc值：0.9945424332338785\n",
      "jpeg.xlsx和不同图片的auc值：0.9965453966789048\n",
      "rotation.xlsx和不同图片的auc值：0.9999995521942556\n",
      "salt_pepper.xlsx和不同图片的auc值：0.9989567269493231\n",
      "speckle.xlsx和不同图片的auc值：0.9945424332338785\n",
      "watermark.xlsx和不同图片的auc值：0.9999996341331792\n",
      "所有相同图片和不同图片 的auc值：0.9984437686265519\n"
     ]
    }
   ],
   "source": [
    "dir_dis_similar_dir(\"../new_data/similar_image_4/\",15,\"./新数据集结果/纯四元数版本_环形分区结果\")\n",
    "dis_different_dir(\"../new_data/different_image/\",15,\"./新数据集结果/纯四元数版本_环形分区结果\")\n",
    "stat_analysis(\"./新数据集结果/纯四元数版本_环形分区结果/\")\n",
    "auc_analysis(\"./新数据集结果/纯四元数版本_环形分区结果/\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "image_qua",
   "language": "python",
   "name": "image_qua"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
